Projeto Final
MAC0460 - Introdução ao aprendizado de máquina
Profª Nina Hirata
A motivação original de nosso projeto é o entendimento dos padrões
climáticos e suas alterações no clima de uma das cidades mais relevantes
num cenário global: Londres. Para isso decidimos fazer um projeto do 1º
formato (utilizando um dataset público), onde a base escolhida
(london_weather.csv, obtida no website Kaggle
através desse
link) pode ser interpretada como uma série temporal, a medida que
mede as variações do tempo (entre 1979 e 2021) de variáveis fixadas em
um único local.
O objetivo do projeto é o uso de um algoritmo de machine learning para que o mesmo consiga prever a categoria de duas das variáveis da database (precipitação e neve) para uma nova informação adicionada ao modelo. As categorias seriam então “chove” ou “não chove” e “neva” ou “não neva”. A ideia será fazer 3 métodos distintos (Regressão Logística, KNN e Árvore de Decisão), compará-los através de suas métricas relacionadas à acurácia e ver qual tem a melhor performance/predição.
A base de dados escolhida foi criada a partir da união de medições oriundas de pedidos de atributos individuais do clima providos pela European Climate Assessment (ECA). As medidas desta base de dados em particular foram gravadas pela estação climática nas redondezas do Aeroporto Heathrow em Londres, Reino Unido. O tamanho original da base de dados escolhida, assim como uma lista dos atributos e suas descrições, está descrito abaixo:
london_weather.csv - 15341 observações x 10
atributos:
date - data em que ocorreu a medição -
(int)
cloud_cover - medição da nebulosidade em oktas -
(float)
sunshine - medição da luz solar em horas (hrs) -
(float)
global_radiation - irradiação medida Watt por metro
quadrado (W/m2) - (float)
max_temp - temperatura máxima registrada em graus
Celsius (°C) - (float)
mean_temp - temperatura média registrada em graus
Celsius (°C) - (float)
min_temp - temperatura mínima registrada em graus
Celsius (°C) - (float)
precipitation - precipitação medida em milímetros
(mm) - (float)
pressure - pressão medida em Pascals (Pa) -
(float)
snow_depth - profundidade da neve medida em
centímetros (cm) - (float)
import pandas as pd
import numpy as np
url = 'https://raw.githubusercontent.com/bmorbin/ML_Project/main/data_raw.csv'
data = pd.read_csv(url)
# Print raw data
print(data.head(5))
# Print count of NA for each column
date cloud_cover sunshine ... precipitation pressure snow_depth
0 19790101 2.0 7.0 ... 0.4 101900.0 9.0
1 19790102 6.0 1.7 ... 0.0 102530.0 8.0
2 19790103 5.0 0.0 ... 0.0 102050.0 4.0
3 19790104 8.0 0.0 ... 0.0 100840.0 2.0
4 19790105 6.0 2.0 ... 0.0 102250.0 1.0
[5 rows x 10 columns]
print(data.isna().sum())
# Removing rows with NA, except the column snow_depth
date 0
cloud_cover 19
sunshine 0
global_radiation 19
max_temp 6
mean_temp 36
min_temp 2
precipitation 6
pressure 4
snow_depth 1441
dtype: int64
data = data.dropna(subset=[col for col in data.columns if col!='snow_depth'])
# Print count of NA again after delete some NA
print(data.isna().sum())
# Applying simple regression model to predict missing values from snow_depth
date 0
cloud_cover 0
sunshine 0
global_radiation 0
max_temp 0
mean_temp 0
min_temp 0
precipitation 0
pressure 0
snow_depth 1418
dtype: int64
from sklearn.linear_model import LinearRegression
def reg_simple(data):
'''
A function to predict snow depth by mean temperature of the day
Returns the maximum snow depth value predict
'''
data_reg_train = data.dropna()
X_reg_train = np.array(data_reg_train.mean_temp).reshape(-1,1)
y_reg_train = np.array(data_reg_train.snow_depth).reshape(-1,1)
reg = LinearRegression()
reg.fit(X_reg_train, y_reg_train)
data_NA = data[data['snow_depth'].isna()]
y_pred = reg.predict(np.array(data_NA.mean_temp).reshape(-1,1))
return(round(max(y_pred)[0],1))
print(f'Minimum snow depth predicted is {reg_simple(data)}')
# So the conclusion is to input 0s in missing values from column snow_depth
Minimum snow depth predicted is 0.2
data.snow_depth[data['snow_depth'].isna()] = 0
# Print count of NA again after input NA values from snow_depth
print(data.isna().sum())
# Convert column date to pandas date type
date 0
cloud_cover 0
sunshine 0
global_radiation 0
max_temp 0
mean_temp 0
min_temp 0
precipitation 0
pressure 0
snow_depth 0
dtype: int64
data['date'] = pd.to_datetime(data['date'].astype(str), format = "%Y%m%d")
# Creating column month
data['month'] = data['date'].dt.month.astype(str)
# Creating column rain indicator
data['rain'] = ["1" if p > 0 else "0" for p in data['precipitation']]
# Creating column snow indicator
data['snow'] = ["1" if p > 0 else "0" for p in data['snow_depth']]
# Setting the target classifier
data['rain_snow'] = data['rain']+data['snow']
# Set column date to index
data = data.set_index('date')
# Print proportion of time wich rains and snows
print(round(data[['snow','rain']].value_counts(normalize=True),3))
snow rain
0 0 0.515
1 0.476
1 1 0.005
0 0.004
dtype: float64
print(round(data.rain_snow.value_counts(normalize=True),3))
# Describe dataframe
00 0.515
10 0.476
11 0.005
01 0.004
Name: rain_snow, dtype: float64
print(data.describe())
cloud_cover sunshine ... pressure snow_depth
count 15261.000000 15261.000000 ... 15261.000000 15261.000000
mean 5.271018 4.345633 ... 101535.436079 0.034336
std 2.067743 4.019266 ... 1049.536452 0.519855
min 0.000000 0.000000 ... 95960.000000 0.000000
25% 4.000000 0.500000 ... 100920.000000 0.000000
50% 6.000000 3.500000 ... 101620.000000 0.000000
75% 7.000000 7.200000 ... 102240.000000 0.000000
max 9.000000 16.000000 ... 104820.000000 22.000000
[8 rows x 9 columns]
Primeiramente devemos encontrar uma solução baseline para o problema proposto. Usando a intuição e o senso comum imagina-se que, em dias muito frios, haverá uma maior probabilidade de nevar, enquanto que dias mais quentes terão maior incidência de chuva e sem neve. Tomemos então isso como nossa solução baseline que será comparada aos modelos categorizados citados acima (Regressão Logística, Árvore de Decisão e KNN).
# 20% for test data
n_test = int(len(data)*0.2)
test_data = data[-n_test:]
# 80% for train and validation
n_train_val = len(data) - int(len(data)*0.2)
df = data[:n_train_val]
target = 'rain_snow'
input_columns = ['cloud_cover','sunshine','global_radiation','max_temp','mean_temp','min_temp','pressure','month']
X, y = df[input_columns], df[target]
# checking if all classes are present in the training set
print(y.value_counts())
00 6314
10 5764
11 67
01 64
Name: rain_snow, dtype: int64
Nos algoritmos a seguir, serão utilizados 30% do conjunto
df (conjunto sem a parte de teste) para validação e 70%
para treino.
Para o algoritmo de classificação por Regressão Logística, serão comparados as diferentes combinações de variáveis.
from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt
import time
from sklearn.linear_model import LogisticRegression
import itertools
start = time.time()
# creating all combinations of the inputs variables
comb = []
for L in range(1, len(input_columns)+1):
for subset in itertools.combinations(input_columns, L):
comb.append(list(subset))
print("Number of combinations: ",len(comb))
Number of combinations: 255
l_total_train = []
l_total_val = []
b=5 #iterations for bootstrap
for c in comb:
l_a_train = []
l_a_val = []
for i in range(1,b):
X_train, X_val, y_train, y_val = train_test_split(X,y,test_size=0.3, random_state=i)
X_train_comb, X_val_comb = X_train[c], X_val[c]
clf = LogisticRegression(multi_class='ovr',solver='lbfgs')
clf = clf.fit(X_train_comb, y_train)
l_a_train += [sum(clf.predict(X_train_comb) == y_train)]
l_a_val += [sum(clf.predict(X_val_comb) == y_val)]
l_total_train += [np.array(l_a_train)]
l_total_val += [np.array(l_a_val)]
mean_train = np.transpose(1-np.array(l_total_train)/len(y_train)).mean(0)
sd_train = np.transpose(1-np.array(l_total_train)/len(y_train)).std(0)
mean_val = np.transpose(1-np.array(l_total_val)/len(y_val)).mean(0)
sd_val = np.transpose(1-np.array(l_total_val)/len(y_val)).std(0)
axis_x = list(range(len(comb)))
# reorder the values by the minimum mean error of train to maximum mean error of train
sort_index = [tup[0] for tup in sorted(list(zip(axis_x,mean_train)),key=lambda tup: tup[1], reverse=True)]
mean_train = mean_train[sort_index]
sd_train = sd_train[sort_index]
mean_val = mean_val[sort_index]
sd_val = sd_val[sort_index]
# index of best combination
best_comb_ID = list(mean_val).index(min(mean_val))
# Plotting the combinations with errors
plt.title('Regressão Logística\nErro x Combinação')
plt.xlabel('ID Combinação')
plt.ylabel('Erro')
plt.plot(axis_x, mean_train, label="Treino")
plt.plot(axis_x, mean_val, label="Validação")
plt.xticks(axis_x, rotation = 90)
plt.fill_between(axis_x, mean_train+sd_train, mean_train-sd_train, alpha=0.5)
plt.fill_between(axis_x, mean_val+sd_val, mean_val-sd_val, alpha=0.5)
plt.axvline(x=best_comb_ID,color='silver', linestyle='dashed')
plt.legend(loc="lower left")
plt.tick_params(axis='x', which='major', labelsize=1)
plt.show()
end = time.time()
print("Spending time to Logistic Regression analysis of combinations:",(end-start)/60, "minutes")
Spending time to Logistic Regression analysis of combinations: 3.7618120551109313 minutes
Para o algoritmo de classificação por Árvore de Decisão, comparou-se os valores de profundidade máxima da árvore que apresentaram menor erro de validação.
from sklearn import tree
l_total_train = []
l_total_val = []
start = time.time()
b=20 # iterator bootstrap (resampling the set)
c=20 # number of complexity based on max_width
for j in range(0,b):
X_train, X_val, y_train, y_val = train_test_split(X,y,test_size=0.3, random_state = j)
l_a_train = []
l_a_val = []
for i in range(1,c):
clf = tree.DecisionTreeClassifier(max_depth = i, criterion = "gini")
clf = clf.fit(X_train, y_train)
l_a_train += [sum(clf.predict(X_train) == y_train)]
l_a_val += [sum(clf.predict(X_val) == y_val)]
l_total_train += [np.array(l_a_train)]
l_total_val += [np.array(l_a_val)]
mean_train = np.transpose(1-np.array(l_total_train)/len(y_train)).mean(1)
sd_train = np.transpose(1-np.array(l_total_train)/len(y_train)).std(1)
mean_val = np.transpose(1-np.array(l_total_val)/len(y_val)).mean(1)
sd_val = np.transpose(1-np.array(l_total_val)/len(y_val)).std(1)
# best max_depth argument to don't overfit
best_depth = list(mean_val).index(min(mean_val))+1
axis_x = list(range(1,c))
# Plotting
plt.title('Árvore de Decisão\nErro x Complexidade')
plt.xlabel('Complexidade: Profundidade Máxima')
plt.ylabel('Erro')
plt.plot(axis_x, mean_train, label="Treino")
plt.plot(axis_x, mean_val, label="Validação")
plt.xticks(axis_x)
plt.fill_between(axis_x, mean_train+sd_train, mean_train-sd_train, alpha=0.5)
plt.fill_between(axis_x, mean_val+sd_val, mean_val-sd_val, alpha=0.5)
plt.axvline(x=best_depth,color='gray', linestyle='dashed')
plt.legend(loc="lower left")
plt.show()
end = time.time()
# time spending to load all algorithm analysis
print("Spending time for Decision Tree analysis of complexity:",(end-start)/60, "minutes")
Spending time for Decision Tree analysis of complexity: 0.43555124203364054 minutes
from sklearn.neighbors import KNeighborsClassifier
b=20 # iterator bootstrap (resampling the set)
c=20 # number of complexity based on number of neighbors
for j in range(0,b):
X_train, X_val, y_train, y_val = train_test_split(X,y,test_size=0.3, random_state = j)
l_a_train = []
l_a_val = []
for i in range(1,c):
clf = KNeighborsClassifier(n_neighbors=c)
clf = clf.fit(X_train, y_train)
l_a_train += [sum(clf.predict(X_train) == y_train)]
l_a_val += [sum(clf.predict(X_val) == y_val)]
l_total_train += [np.array(l_a_train)]
l_total_val += [np.array(l_a_val)]
mean_train = np.transpose(1-np.array(l_total_train)/len(y_train)).mean(1)
sd_train = np.transpose(1-np.array(l_total_train)/len(y_train)).std(1)
mean_val = np.transpose(1-np.array(l_total_val)/len(y_val)).mean(1)
sd_val = np.transpose(1-np.array(l_total_val)/len(y_val)).std(1)
best_nn = list(mean_val).index(min(mean_val))+1
axis_x = list(range(1,c))
# Plotting for view the best number of neighbors
plt.title('KNN\nErro x Complexidade')
plt.xlabel('Complexidade: número de vizinhos')
plt.ylabel('Erro')
plt.plot(axis_x, mean_train, label="Treino")
plt.plot(axis_x, mean_val, label="Validação")
plt.xticks(axis_x)
plt.fill_between(axis_x, mean_train+sd_train, mean_train-sd_train, alpha=0.5)
plt.fill_between(axis_x, mean_val+sd_val, mean_val-sd_val, alpha=0.5)
plt.axvline(x=best_nn,color='silver', linestyle='dashed')
plt.legend(loc="lower left")
plt.show()
end = time.time()
print("Spending time for KNN algorithm analysis:",(end-start)/60, "minutes")
Spending time for KNN algorithm analysis: 3.720773720741272 minutes
A partir da seleção de parâmetros de cada algoritmo realizada nas seções anteriores, compara-se então os modelos entre si para analisar qual possui melhor acurácia.
from sklearn.metrics import plot_confusion_matrix
X_train, X_val, y_train, y_val = train_test_split(X,y,test_size=0.3, random_state = 0)
X_train_comb, X_val_comb = X_train[comb[best_comb_ID]], X_val[comb[best_comb_ID]]
print("Accuracy of Regression Model: ",sum(LogisticRegression(multi_class='ovr',solver='lbfgs').fit(X_train_comb, y_train).predict(X_val_comb) == y_val)/len(y_val))
Accuracy of Regression Model: 0.6033306033306033
plot_confusion_matrix(LogisticRegression(multi_class='ovr',solver='lbfgs').fit(X_train_comb, y_train), X_val_comb, y_val)
<sklearn.metrics._plot.confusion_matrix.ConfusionMatrixDisplay object at 0x0000000049EB7EB0>
plt.show()
print("Accuracy of Decision Tree: ",sum(tree.DecisionTreeClassifier(max_depth=best_depth, criterion = "gini").fit(X_train, y_train).predict(X_val) == y_val)/len(y_val))
Accuracy of Decision Tree: 0.7447447447447447
plot_confusion_matrix(tree.DecisionTreeClassifier(max_depth=best_depth, criterion = "gini").fit(X_train, y_train), X_val, y_val)
<sklearn.metrics._plot.confusion_matrix.ConfusionMatrixDisplay object at 0x000000004AA3CFD0>
plt.show()
print("Accuracy of KNN: ",sum(KNeighborsClassifier(n_neighbors=best_nn).fit(X_train, y_train).predict(X_val) == y_val)/len(y_val))
Accuracy of KNN: 0.6743106743106743
plot_confusion_matrix(KNeighborsClassifier(n_neighbors=best_nn).fit(X_train, y_train), X_val, y_val)
<sklearn.metrics._plot.confusion_matrix.ConfusionMatrixDisplay object at 0x0000000005CDB280>
plt.show()
Conclui-se então que o algoritmo árvore de decisão apresenta maior acurácia, além de ser o algoritmo mais rápido para análise. Portanto, roda-se a seguir o modelo com o conjunto de dados de teste:
# Fit the model selected
model = tree.DecisionTreeClassifier(max_depth=best_depth, criterion = "gini").fit(X_train, y_train)
X_test = test_data[input_columns]
y_test = test_data[target]
print("Accuracy of Decision Tree to test data: ",sum(model.predict(X_test) == y_test)/len(y_test))
Accuracy of Decision Tree to test data: 0.7260812581913499
plot_confusion_matrix(model, X_test, y_test)
<sklearn.metrics._plot.confusion_matrix.ConfusionMatrixDisplay object at 0x000000004C7362B0>
plt.show()
Todos códigos referentes ao projeto podem ser encontrados nesse repositório.